# load and attach packages----
library(ggplot2)
library(margins)
library(plm)
library(haven)
library(jtools)
library(Zelig) # Can be used for rare events logit, but cannot handle robust clustered errors :(
library(bife) # For fixed effects logit
library(pcse)
library(visreg) # Plots Interactions and predicted values against residuals
library(effects) # Plots effects including interactions
library(tidyverse) # For data manipulation
library(dplyr) # For data manipulation
library(Hmisc) # for summary statistics
library(foreign) # allows the exporting of R back to Stata

# These packages are for making MSWord tables
library(sjPlot)
library(sjmisc)
library(sjlabelled)
library(devtools)
library(PanelMatch)
library("car")
library(tidyr)

# Load Data----
SaCrat_data01.2 <- read_dta("SanctionsAndStockMarkets.dta")

# convert SD and countyrmo to integers. findAllTreated needs them to be integers later
SaCrat_data01.2$ccode <- as.integer(SaCrat_data01.2$ccode)
SaCrat_data01.2$countyrmo <- as.integer(SaCrat_data01.2$countyrmo)
SaCrat_data01.2<-subset(SaCrat_data01.2, (!is.na(SaCrat_data01.2$SG20D1)))
SaCrat_data01.2 <- data.frame(SaCrat_data01.2)


#PanelMatch 
  # With CBPSW refinement 
  PM.resultsG.CBPSW <- PanelMatch(lag = 12, time.id = "countyrmo", unit.id = "ccode", 
                            treatment = "SG20DT5", refinement.method = "CBPS.weight", 
                            data = SaCrat_data01.2, match.missing = T, 
                            covs.formula = ~ I(lag(StMk, 1:3)) + IntCon +CW +prelec +ideoLef +openk +EXRat +CPI +lnTR +lnGDP + 
                              nonSG20D +Dev_Cou, 
                            size.match = 10, qoi = "att" ,outcome.var = "StMk", verbose=TRUE, 
                            lead = 0:3, forbid.treatment.reversal = FALSE,
                            use.diagonal.variance.matrix = TRUE)
  # With no refinement
  PM.resultsG.none <- PanelMatch(lag = 12, time.id = "countyrmo", unit.id = "ccode", 
                            treatment = "SG20DT5", refinement.method = "none", 
                            data = SaCrat_data01.2, match.missing = T, 
                            covs.formula = ~ I(lag(StMk, 1:3)) + IntCon +CW +prelec +ideoLef +openk +EXRat +CPI +lnTR +lnGDP + 
                              nonSG20D +Dev_Cou, 
                            size.match = 10, qoi = "att" ,outcome.var = "StMk", verbose=TRUE, 
                            lead = 0:3, forbid.treatment.reversal = FALSE,
                            use.diagonal.variance.matrix = TRUE)

# generate panel estimates With CBPS.weight Refinement
  PE.results95G.CBPSW <- PanelEstimate(inference = "bootstrap", sets = PM.resultsG.CBPSW, 
                                      data = SaCrat_data01.2, CI = .95, ITER = 1000, set.seed(515))
  summary(PE.results95G.CBPSW)

  PE.results90G.CBPSW <- PanelEstimate(inference = "bootstrap", sets = PM.resultsG.CBPSW, 
                                      data = SaCrat_data01.2, CI = .90, ITER = 1000, set.seed(515))
  summary(PE.results90G.CBPSW)
 
#################################
  # Plot ATT with 90% and 95% standard errors -- CBPSW
#################################
  results95G.CBPSW <- as.data.frame(summary(PE.results95G.CBPSW)$summary)
  results90G.CBPSW <- as.data.frame(summary(PE.results90G.CBPSW)$summary)

theme_set(theme_bw(base_size = 15))

# PLOT FIVE YEAR WINDOW
G <- ggplot(data = results95G.CBPSW, aes(y = estimate, x = 1:4)) +
  geom_errorbar(aes(ymin = `2.5%`, 
                    ymax = `97.5%`), width=0, size = 1) +
  geom_errorbar(data = results90G.CBPSW, aes(ymin = `5%`,
                                      ymax = `95%`), 
                width=0, size = 2) +
  geom_point(size = 6.) +
  geom_hline(yintercept = 0, linetype = "dashed", size = 0.5) +
  scale_y_continuous() +
  scale_x_discrete(name="Months Since G-20 Sanction", limits=c("0"="t+0","1"="t+1","2"="t+2","3"="t+3"))+
  theme(panel.grid.major.x = element_blank(), panel.grid.minor.x = element_blank(),
      panel.grid.minor.y = element_blank()) +
  ylab("Estimated Effect of Treatment on the Treated")

G

#################################
# CHECK AND PLOT CHANGE IN COVARIATE BALANCE
#################################
# Calculate covariate balance
refined_bal <- as.data.frame(get_covariate_balance(PM.resultsG.CBPSW$att, SaCrat_data01.2, 
                                                   covariates = c("IntCon", "CW", "prelec", "ideoLef", "openk",
                                                                  "EXRat", "CPI", "lnTR", "lnGDP", "Dev_Cou", "nonSG20D"), 
                                                   plot = F))
# Reformat for plotting
refined_bal <- gather(refined_bal, key = "covariate", value = "ref_sd")
refined_bal$t <- 12:0
refined_bal$abs_ref_sd <- abs(refined_bal$ref_sd)

# Same procedure for "unrefined" sets
# Calculate covariate balance 
unrefined_bal <- as.data.frame(get_covariate_balance(PM.resultsG.none$att, SaCrat_data01.2, 
                                                     covariates = c("IntCon", "CW", "prelec", "ideoLef", "openk",
                                                                    "EXRat", "CPI", "lnTR", "lnGDP", "Dev_Cou", "nonSG20D"), 
                                                     plot = F))
# Reformat for plotting
unrefined_bal <- gather(unrefined_bal, key = "covariate", value = "unref_sd")
unrefined_bal$t <- 12:0
unrefined_bal$abs_unref_sd <- abs(unrefined_bal$unref_sd)

# Merge into a single dataframe
balance_test <- merge(refined_bal, unrefined_bal, by = c("covariate", "t"))
balance_test$covariate <- as.factor(balance_test$covariate)

# Create a reference line
refLine <- data.frame(x = c(0, 1), y = c(0, 1))

# Plot the covariate balance of "refined" set against "unrefined" set
ggplot(data = balance_test, 
       aes(x = abs_unref_sd, y = abs_ref_sd, color = covariate)) +
  geom_point(size = 4) +
  scale_y_continuous(limits = c(0, 1), breaks = c(0, .25, .5, .75, 1)) +
  scale_x_continuous(limits = c(0, 1), breaks = c(0, .25, .5, .75, 1)) +
  geom_line(data = refLine, aes(x = x, y = y, color = NULL), linetype = 2) +
  ylab("Post-Refinement") +
  xlab("Pre-Refinement")

rm(refLine, refined_bal, unrefined_bal, balance_test)
